// @author Maciej Chałapuk (maciej@chalapuk.pl)
// vim: ts=2 sw=2 expandtab

#ifndef STOCHASTIC_NUMBER_HPP
#define STOCHASTIC_NUMBER_HPP

#include "bitset.hpp"

#include <iostream>
#include <cstdint>
#include <cassert>
#include <random>


namespace experiment {


template <
  size_t bitsize_,
  class bitset_ = bitset<bitsize_>,
  class seedfunc_ = std::random_device
  >
class stochastic_number {
public:
  static const size_t c_bitSize = bitsize_;
  static constexpr long double c_epsilon = (long double) 1 / c_bitSize;

  typedef bitset_ bitset;
  typedef seedfunc_ seedfunc;
  typedef stochastic_number<c_bitSize, bitset, seedfunc> type;

  stochastic_number() {
    m_bitset.reset();
  }
  template <class number>
  stochastic_number(number num) {
    assign(num);
  }
  stochastic_number(stochastic_number const& rhs) {
    assign(rhs.to_double());
  }
  stochastic_number(stochastic_number && rhs)
    : m_bitset(std::move(rhs.m_bitset)) {
  }
  template <class number>
  stochastic_number const& operator=(number num) {
    assign(num);
    return *this;
  }
  stochastic_number const& operator *=(stochastic_number const& rhs) {
    m_bitset &= rhs.m_bitset;
    return *this;
  }
  stochastic_number operator*(stochastic_number const& rhs) const {
    return (clone() *= rhs);
  }
  template <class number>
  stochastic_number const& operator +=(number rhs) {
    assign(to_double() + rhs);
    return *this;
  }
  template <class number>
  stochastic_number operator+(number rhs) const {
    return (clone() += rhs);
  }
  template <class number>
  stochastic_number const& operator /=(number rhs) {
    assign(to_double() / rhs);
    return *this;
  }
  template <class number>
  stochastic_number operator/(number rhs) const {
    return (clone() /= rhs);
  }
  template <class number>
  stochastic_number const& operator -=(number rhs) {
    assign(to_double() - rhs);
    return *this;
  }
  template <class number>
  stochastic_number operator-(number rhs) const {
    return (clone() -= rhs);
  }
  double to_double() const {
    return (m_bitset.count() * c_epsilon);
  }
  operator double() const {
    return to_double();
  }

private:
  bitset m_bitset;

  template <class number>
  void assign(number num);
  size_t dist(size_t a, size_t b) const {
    return static_cast<size_t>(
        std::abs(static_cast<ssize_t>(a) - static_cast<ssize_t>(b))
        );
  }
  size_t min(size_t a, size_t b, size_t c, size_t d) const {
    return std::min(a, std::min(b, std::min(c, d)));
  }
  stochastic_number clone() const {
    stochastic_number retval;
    retval.m_bitset = m_bitset;
    return retval;
  }

  template <size_t bitSize, class bitset, class seed>
  friend std::ostream& operator<<(std::ostream &,
      stochastic_number<bitSize, bitset, seed> const&);
};

template <size_t bitsize, class bitset, class seedfunc>
template <class number>
void stochastic_number<bitsize, bitset, seedfunc>::assign(number num) {
  assert(num >= 0);
  assert(num <= 1);

  size_t curBitCount = m_bitset.count();
  size_t numBitCount = static_cast<size_t>(num / c_epsilon);
  size_t flipBitCount = c_bitSize - curBitCount;

  const size_t distance = dist(curBitCount, numBitCount);
  const size_t flipDistance = dist(flipBitCount, numBitCount);
  const size_t beginDistance = numBitCount;
  const size_t endDistance = dist(c_bitSize, numBitCount);

  size_t changeCount = 0;
  size_t minDist = min(distance, flipDistance, beginDistance, endDistance);
  if (minDist == distance) {
    changeCount = distance;

  } else if (minDist == flipDistance) {
    changeCount = flipDistance;
    curBitCount = flipBitCount;
    m_bitset.flip();

  } else if (minDist == beginDistance) {
    changeCount = beginDistance;
    curBitCount = 0;
    m_bitset.reset();

  } else if (minDist == endDistance) {
    changeCount = endDistance;
    curBitCount = c_bitSize;
    m_bitset.reset();
    m_bitset.flip();
  }

  static seedfunc seed;
  static std::mt19937 generator(seed());
  static std::uniform_int_distribution<> distribution(0, c_bitSize-1);

  size_t i = 0;
  bool val = (curBitCount < numBitCount);

  while (i != changeCount) {
    size_t idx = distribution(generator);
    if (m_bitset[idx] != val) {
      m_bitset[idx] = val;
      i += 1;
    }
  }
}

template <size_t bitsize, class bitset, class seedfunc>
std::ostream& operator<<(std::ostream &out,
    stochastic_number<bitsize, bitset, seedfunc> const& num) {
  return out << static_cast<double>(num);
}


typedef stochastic_number<1024> stochastic1024_t;
typedef stochastic_number<2048> stochastic2048_t;
typedef stochastic_number<4096> stochastic4096_t;
typedef stochastic_number<8192> stochastic8192_t;
typedef stochastic_number<16384> stochastic16384_t;
typedef stochastic_number<32768> stochastic32768_t;
typedef stochastic_number<65536> stochastic65536_t;
typedef stochastic_number<131072> stochastic131072_t;
typedef stochastic_number<262144> stochastic262144_t;
typedef stochastic_number<524288> stochastic524288_t;
//typedef stochastic_number<> stochastic_t;

} // namespace experiment


#endif // include guard

